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Motivated by a recent work on the metabolism of carbohydrates in bacteria, we study the kinetics and thermo¬ 
dynamics of two classic models for reversible polymerization, one preserving the total polymer concentration 
and the other one not. The chemical kinetics is described by rate equations following the mass-action law. We 
consider a closed system and nonequilibrium initial conditions and show that the system dynamically evolves 
towards equilibrium where detailed balance is satisfied. The entropy production during this process can be ex¬ 
pressed as the time derivative of a Lyapunov function. When the solvent is not included in the description and the 
dynamics conserves the total concentration of polymer, the Lyapunov function can be expressed as a Kullback- 
Leibler divergence between the nonequilibrium and the equilibrium polymer length distribution. The same result 
holds true when the solvent is explicitly included in the description and the solution is assumed dilute, whether 
or not the total polymer concentration is conserved. Furthermore, in this case a consistent nonequilibrium ther¬ 
modynamic formulation can be established and the out-of-equilihrium thermodynamic enthalpy, entropy and 
free energy can be identified. Such a framework is useful in complementing standard kinetics studies with the 
dynamical evolution of thermodynamic quantities during polymerization. 

PACS numbers: 05.70.Ln, 05.40.-a, 05.70.-a 


I. INTRODUCTION 


The processes of aggregation or polymerization are ubiquitous in Nature, for instance they are present in the poly¬ 
merization of proteins, the coagulation of blood, or even in the formation of stars. They are often modeled using the 
classic coagulation equation derived by Smoluchowsky [1,2]. In the forties, P. Flory [3,4] developed his own approach 
for reactive polymers, with an emphasis on their thermodynamic properties and on their most likely (equilibrium) size 
distribution. At this time, only irreversible polymerization was considered. The first study of the kinetics of reversible 
polymerization combining aggregation and fragmentation processes was carried out by Tobolsky et al. [5]. 

In the seventies, reversible polymerization became a central topic in studies on the association of amino acids into 
peptides and on the self-assembly of actin. The thermodynamics of assembly of these polymers forms the topic of the 
now classical treaty by F. Oosawa and S. Asakura [6]. At about the same time, T. L. Flill made many groundbreaking 
contributions to non-equilibrium statistical physics and thermodynamics, which allowed him to describe not only the 
self-assembly of biopolymers like actin and microtubules, but also to address much more complex questions such that 
of free energy transduction by biopolymers and complex chemical networks [7]. In the eighties, R. J. Cohen and G. 
B. Benedek [8] revisited the work of Flory and Stockmayer, by showing that the Flory polymer length distribution 
is obtained under the assumption of equal free energies of bond formation for all bonds of the same type. They also 
showed that the kinetically evolving polymer distribution does not have the Flory form in general, and they analyzed 
the irreversible kinetics of the sol-gel transition. 

More recently, the specific conditions on the kernels of aggregation and fragmentation, for which equilibrium solu¬ 
tions of the Hory type exist have been analyzed [9-11]. When these conditions are not met, reversible polymerization 
models admit interesting nonequilibrium phase transitions which are beginning to be investigated [12, 13]. 

The kinetic rate equations of reversible polymerization have broad applications. For instance, in studies on the 
origin of life, these equations describe the appearance of long polymer chains in the primordial soup [14]. These 
equations are also used to describe the formation of protein clusters in membranes [15-17] or the self-assembly of 
carbohydrates (also called glycans) [18]. In the latter case, a very large repertoire of polymer structures and enzymes 
are involved in the synthesis and degradation of these polymers. Since it is hardly possible to model all the involved 
chemical reactions, the authors of this work, Kartal et al., introduced a statistical approach to explain experiments 
which they have performed using mixtures of such polymers with the appropriate enzymes. Their study underlines the 
importance of entropy as a driving force in the dynamics of these polymers: under its action a monodisperse solution 
of such biopolymers, which is placed in a closed reactor with the appropriate enzymes, typically admits an exponential 
distribution of polymer length as equilibrium distribution, in agreement with maximum entropy arguments used by R 
Flory [3,4]. 

This recent work of Kartal et al. [18], motivated us to construct appropriate dynamics which converge on long 
times towards such equilibrium distributions. In order to complement this with an analysis of the time evolution of 
thermodynamic quantities, we rely on stochastic thermodynamics (for general reviews see [19-22]). While this recent 
branch of thermodynamics has been used extensively in the literature for chemical reaction networks [23-30] and 
copolymerization processes [31] at the level of the stochastic chemical master equation, its application to the level of 
mean-field kinetic rate equations is more recent [32]. 

In this paper, we precisely use this level of description based on mean-field rate equations. Implicitly, we assume 
reaction-limited polymerization. Naturally, if the reactions are too fast or the mobility of the polymers too low, a 
mean-field approach may not be sufficient and diffusion processes should be accounted for. We focus on two main 
models of reversible polymerization which reproduce the equilibrium distributions found in Ref. [18]: In the first one, 
there is only one conservation law (the total number of monomers) while in the second one, there are two (the total 
number of monomers and of polymers). 

The plan of the paper is as follows. In section II, we study reversible polymerization using general rate equations 
compatible with one conservation law (the total number of monomers). This is done first at the one-fluid level for 
which there is no solvent, and then at the two-fluid level, for which there is a solvent. In section III, we apply this 
general framework to a specific model called String model, in which the rates of aggregation and fragmentation are 
constant. In section IV, we extend the previous case of reversible polymerization with a single conservation law to 
the case where there are two conservation laws, namely the total number of monomers and of polymers. In section V, 
inspired by Ref. [18], we study two specific examples of reversible polymerization with two (resp. three) conservation 
laws, namely the kinetics of glucanotransferases DPEl (resp. DPE2). Eor both cases, we construct the dynamics 
which converge towards the equilibrium distributions found in Ref. [18], and we discuss their properties from the 
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standpoint of nonequilibrium thermodynamics. 


II. REVERSIBLE POLYMERIZATION WITH ONE CONSERVATION LAW 


We consider a reversible polymerization process made of the following elementary reactions 

[n] + [m] ^[n + m], (1) 


where the forward and backward reaction rates, namely k~^ and are functions of the polymer lengths n and 
m, which are strictly positive and symmetric under a permutation of n and m. We denote by c; the concentration of 
polymers of length 1. The evolution of this quantity is ruled by 
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( 2 ) 


which preserves the total concentration of monomers (i.e. the concentration that one would get if all the polymers 
were broken into monomers) M = concentration of polymers c = Therefore in 

this case, there is only one conservation law, that of M, and one can assume M = 1. Note also that Eq. 2 generalizes 
the Becker-Doring equations which describe the dynamic evolution of clusters that gain or loose only one unit at a 
time [33, 34] 

Assuming that the reactions (1) can be treated as elementary, the entropy production rate is given by [35, 36] 

S = I In > 0, (3) 

^ KnmCn+m 


where R is the gas constant. The entropy production rate vanishes when the system reaches equilibrium, i.e. when 
and only when detailed balance is satisfied; 
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Using the inequality In a; < a; — 1, which holds for all a; > 0, one easily proves that the following quantity is 
non-negative, convex and vanishes only at equilibrium 

L = c/In- i?(c-> 0. (5) 


Indeed, by taking the time derivative of L{t) and using the definition of c, one obtains dL/dt = RJ2i logcj/c^®. 
Now using (2), we obtain two terms. The first term is 
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By adding these two contributions and using the detailed balance condition of (4) as well as (3), one finds that 

E = -i?^Qln^ = -^L>0. (6) 

i 

This proves that L is a Lyapunov function, i.e. a non-negative, monotonically decreasing function, which vanishes 
at and only at equilibrium. The existence of such a function implies that the dynamics will always relax to a unique 
equilibrium state. In [37], the authors show that for a chemical system containing a finite number of homogeneous 
phases, a Gibbs free energy function exists that is minimum at equilibrium. In the context of reacting polymers, a 
similar free energy function has been derived in [9]. In view of the non-increasing property of this function, the 
authors have coined the term “F-Theorem”. We discuss below its relation to the more usual H theorem. 
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A. One-fluid version 


Until now, our model was exclusively expressed in terms of concentrations and could be used to study non-chemical 
dynamics such as population dynamics. We now introduce the one fluid model, where the solvent is not described 
either by choice or because it is absent. The molar fractions of the polymers of length I, with I > 1, are xi = c;/c 
while the Lyapunov function is 


L = + Rc\n— -i?(c — > 0. (7) 

/>! ^*=9 

It is important to note that this Lyapunov function is in general distinct from the relative entropy or Kullback-Leibler 
divergence between the distribution xi and which represents only the first term in Eq. (5). The reason for this 
difference is that c;, contrary to xi, cannot always be interpreted as a probability since its norm is not always conserved. 
The two quantities become however equivalent (i.e. they only differ by a constant) when the total concentration of 
polymers c{t) is constant in time. If furthermore xi is constant in time (the meaning of this assumption will 
become clear in the two fluid model), then the negative of the Shannon entropy (which is related to the classic notion 
of free energy of mixing introduced in [3, 4] as explained in [38]) 

Ssh = -R ^ xi 111 xi , ( 8 ) 

^>1 

becomes a Lyapunov function and (6) reduces to the famous Boltzmann “H theorem”. 


B. Two-fluid version 

In order to make contact with thermodynamics, we now introduce the two fluid model which includes the solvent ex¬ 
plicitly in the list of chemical species. Then, the molar fraction of the polymer of length / becomes yi (t) = ci (t) jCit'), 
which importantly is now defined with respect to the total concentration of all species including the (time-independent) 
solvent concentration cq (water for instance): C{t) = J2i>o = c(t) -I- cq, where c{t) = J2i>i ci{t). If there is 
no solvent, cq = 0 and one recovers the molar fraction of the previous section which was denoted by xi. In dilute 
solution, since yo is very close to one and the other yi are much smaller, C{t) becomes almost constant: C « cq. The 
chemical potential of a polymer of length I in a dilute solution is defined by yi = + RT In yi, where — Ts^ 

is the standard reference chemical potential and and s° are the standard enthalpy and entropy respectively. We 
restrict ourselves here to ideal solutions, and by this we assume that this form of chemical potential applies not only 
to the polymers (I ^ 0) but also to the solvent. An interesting study of the effect of non-ideality on the time evolution 
thermodynamic quantities during reversible polymerization can be found in Ref. [39]. 

Let us define the intensive enthalpy function as 


H = J2yihl 

l>0 


(9) 


and the entropy function of this two-fluid model as 

5 = ^yi{s^ - Rlnyi). (10) 

/>o 

Their extensive counterparts are TL = CH and S = CS. In Eq. (10), the first term proportional to s° therefore 
represents the entropic contribution due to the disorder in the internal degrees of freedom of each polymer, while the 
second term represents the nonequilibrium entropy in the distribution of the variables yi [26]. Let us introduce the 
intensive free enthalpy 


G = H-TS = Y,yiyi, 

l>0 


( 11 ) 
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where we used the definition of the chemical potential in the last equality, and its extensive counterpart G = CG. 

Since the change of chemical potential associated to each reaction must vanish at equilibrium, i.e. = ^n+m — 
Mn ~ Mm = 0, using the definition of the chemical potential, we find that 
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Combining (12) with (4), we obtain that the kinetic constants must satisfy local detailed balance 


( 12 ) 
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-(M„+m - M„ - Mm): 


(13) 


where denotes the equilibrium value taken by C. We note that since the chemical reactions do not involve the 
solvent, we formally define the rate constants with any zero subscript (n or to) to be zero. Naturally, Eq. (13) is not 
applicable in this case. The validity of Eq (13) relies on two main assumptions: the first one is that of dilute solutions, 
while the second one is the ideality of the heat bath. The latter assumption means that there are no hidden degrees 
of freedom which can dissipate energy in the chemical reactions under consideration, which implies in particular that 
these reactions must be elementary. 

Using the definitions of the enthalpy and entropy functions, we find that 
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As a result, one finds that 
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where we used Eq. (12) to obtain the last equality. By including the solvent in the sum, as the term corresponding to 
Z = 0, we have J2i>o yi — ^ therefore J2i>o Vi- ~ rewrite the above equation as 


dG 


dt 
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(17) 


The entropy production (3), using (13) and the chemical potential definition, may be written as 

2^ _^ (j 

TE = —— ^ ~ k^^Cn+m) {d-n+m ~ Mn ~ Mm ~ RT \ll ) — 0 * 

n,m 


(18) 


Using Eq. (14) and Eq. (15), it can be rewritten as 


TS =- 

dt 


n-TS- RTC In ^ - C) 


> 0 . 


(19) 


The first term is the heat flow, the second the entropy change and the third and fourth terms represent a contribution 
due to the change in the total concentration. It is important to note that within the two-fluid model with ideal solutions, 
since C{t) is essentially constant, these latter two contributions are negligible. Neglecting these terms, the entropy 
production can be expressed as the change in free energy which is also equal to a change in Kulback-Leibler divergence 
between the nonequilibrium and the equilibrium polymer distribution 


rs = - 


dt 



( 20 ) 
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We finally note that when all the polymerization reactions are neutral from a standard chemical potential standpoint 
i.e. when fJ-n+m = + Mm has that 

t = >”!«). pi) 

/>o i>0 

which together with Eq. (17) implies dt{C Ui = 0- Since C can be assumed constant, the Lyapunov function 

can be expressed only in terms of the Shannon entropy constructed from yi instead of the full KL divergence. The 
dynamics can then be compared to a Boltzmann equation where the relaxation to equilibrium is purely driven by the 
maximization of the Shannon entropy, as in the H-Theorem. 


III. APPLICATION TO THE STRING MODEL 


As a simple realisation of the reversible polymerization given by Eq. (1), we now consider the String model 
which assumes constant rates of aggregation and fragmentation, independent of the length of the reacting polymers. 
Eollowing [13], we choose = 2 and k~^ = 2A. Erom Eq. (2), the dynamics follows 

Cl = CiCj + Cj — 2c/c — X{1 — l)ci, (22) 

i+j=l j>l 

where we used the fact that (Z —1) combinations of z and j satisfy the relation z+j = 1. The detailed balance condition 
defining equilibrium implies that = Ac®^^ , which admits one parameter solutions of the form = X/3K 

Assuming the total monomer concentration to be M = = 1, one finds 

/3 = 1 + - - Y^A + —, (23) 

since the solution ci must decay at large 1. 


A. One-fluid version of the String model 

The evolution of the length distribution in the String model can be obtained as a function of time by explicit numer¬ 
ical integration. The results, starting from time f = 0, are shown in Eig. 1. 



FIG. 1. Evolution of the size distribution of polymers c; for the String model at different times t = 0, t = 0.2, t = 0.5 and 
t = 1 starting with a monodisperse distribution with characteristic length Z = 4 at time f = 0. The transition rates are constant and 
correspond to A = 1. A rapid convergence towards an exponential distribution can already be seen after a time f ~ 1. 


The evolution equation (22) can be solved using an exponential ansatz of the form [13] 

c/(f) = (1 - a(f))^a'"^(f), (24) 
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FIG. 2. Left: Shannon entropy of the one fluid version of the String model, Ssh, as a function of time. The various curves represent 
different initial conditions of the form ci{t = 0) = Sim, for different total monomer concentration M = 1, 2, 3 and 4. Right: Total 
polymer concentration c{t) for the same initial conditions, as found by numerical integration of kinetic equations (solid line) or 
from an exact expression using Eq. (36) (symbols). Both figures have been done with the same rates as Fig. 1. 


which satisfies the conservation of total number of monomers. The resulting differential equation for a(t) is d = 
(1 — a)^ — Xa. This equation can be easily solved with a monomer-only initial condition, ci = Si^i, which translates 
into the condition a(0) = 0. Unfortunately, the exponential ansatz cannot be used to describe more general initial 
conditions, which cannot be accounted for by such a simple (-independent condition on a(0). For the monomer-only 
initial condition, the following explicit solution is obtained: 

l-exp{-{a+-a-)t} 

a[t) = - - — - -^ 

Q+ — a_ exp { —(a+ — a_)t} 

a± = l + -±yA+—, 

where the two roots a_ and a+ are related by a-a+ = 1. At long times, the RHS of Eq. (25) tends towards 
l/a+ = a- = P, so that the system approaches the equilibrium distribution c/(oo) = = (1 — 

Using (25), one can obtain the explicit time evolution of the quantities of interest; The total polymer concentration 
is time-dependent and reads 


(25) 

(26) 


c(f) = ^ C; (f) = 1 - a(t). (27) 

/>1 

The Shannon entropy at the one fluid level, and defined in Eq. (8), is given by 

=-ln(l - a(f)) - Ina(f), (28) 

and reaches its equilibrium value Ssh{oo) = — In (1 — — (a_/(1 — a-)) In a_ for long times. Its rate of change 

is 


Ssh = ~o(f) Ina(f) -F A 


a{t) 




■ Ina(f), 


(29) 


while the entropy production rate given by Eq. (6) is 

S = i? ("2 In - In ((1 - a(f))2 - Aa(f)) . 

\ 1 — q;_ a- J 


(30) 


Eor completeness, we also discuss an approach using generating functions to study the String model without resort¬ 
ing to the exponential ansatz which is restricted to the monomer-only initial condition. The full dynamics remains 
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nevertheless complex to solve within this approach. Introducing the generating function t) 
using Eq. (22), we obtain the dynamical equation 


dt 


= X^ -2cx + 2Xj2Y.CjZ^ - \iz^-x) , 


i j>i 


which can be simplified since for \z\ < 1, ^j>i ~ ~ x)/(l ~ so that 


ot \ — z 




Y.i>iCi{t)z^ and 

(31) 

(32) 


The stationary solution of this differential equation, i.e. the solution of dx/dt = 0 has the following form: 

2 

X^'^iz) = x(.Z,t ^ oo) = ---r-(33) 

which satisfies x(l,f —>■ oo) = and oo) = c®'^(A + c®'^)/A. Since M = = 

X^(l, t —>■ 00 ) = 1, one recovers the equilibrium state obtained before, namely = 1 — /3, with (3 given by Eq. (23). 
Eor an arbitrary M, one deduces from Eq. (33) that the equilibrium polymer length distribution is = Aa^ where 
a = I (A + d'^). With our choice of initial condition of the form c/(0) = &im, we have J2i>i Therefore, 

the equilibrium solution is 


c 


eq 

I 


{d'if 

M 



(34) 


Besides the stationary solution, it is also possible to obtain analytically the evolution of the total concentration c(f). 
To show this, we take the limit 2 —1 in Eq. (32). Using I’Hospital rule, we obtain 

c = -c^ + A(M - c). (35) 

The explicit solution of Eq. (35) for c(0) 

5 im is 

c{t) = 


= 1, as imposed by our choice of initial conditions of the form ci{t = 0) = 


A , A + 2\ a 

— tanh - A + arctanh —^—-. 

2 V2 A I 2 


(36) 


where A = A(A + AM). One can verify that Eq. (27) is recovered for the monomer only initial condition, M = 1, 
as expected. Eurthermore, one recovers that c{t) tends towards = (A — A)/2 as t —oo, which is the equilibrium 
concentration entering in Eqs. (33) and (34). 

Incidentally, one may wonder whether this behavior of the total concentration agrees with the predictions of Ref. 
[6] regarding the notion of critical concentration in reversible polymerization. This is indeed the case: if one evaluates 
the concentration of monomers ci as a function of M, with Eq (34) and eliminating c®® using the above expression, 
one finds a function of M which first increases rapidly and then reaches a plateau for M > X. Naturally, this is not a 
sharp transition but rather a cross-over between two regimes. One can also look at the average length of the polymer 
M / d‘^ which increases significantly when M becomes larger than A. Both features indicate that A represents the 
critical concentration of this model [6]. 

As shown in the right part of Eig. 2, which has been obtained by explicit numerical integration, the polymer 
concentration c{t) either decreases as a function of time for the monomer only initial condition (M = 1) or increases 
as a function of time when the initial condition corresponds to polymers of length 3 or above (M > 3), while it remains 
constant for the case of dimers (M = 2). Intuitively, when the initial condition is monomer-only (M = 1), there is 
mainly aggregation of monomers, so that the net concentration must decrease with time. On the other hand, if the 
initial solution consists of long polymers (M > 2), the probability of fragmentation is higher than that of aggregation. 
As a result, the total concentration must increase with time. Eor dimer-only initial condition (M = 2), the probabilities 
of fragmentation is same as that of aggregation, so that the net concentration stays constant. As a result of this time 
dependence of c(f), we also see in the left part of Eig. 2, that the Shannon entropy Ssh does not always increase 
monotonically as a function of time. It does so for M <2 but not for M > 2, where it presents an overshoot before 
reaching its equilibrium value. Such an overshoot reveals that the Shannon entropy is not a Lyapunov function L as 
discussed in the previous section. 
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B. Two-fluid version of the String model 


One of the main difference between the two fluid approach as compared with the previous one with a single fluid, 
is the existence of the local detailed balance condition, namely Eq. (13), which connects the rate constants to the 
difference of standard chemical potentials. Further, the specific form of the standard chemical potentials enters in the 
equilibrium length distribution of the polymers and in the kinetics of the self-assembly process. 

For instance, if the polymers self-assemble linearly, the standard chemical potential of a polymer of length I, for 
I > 1, may be written as = — (Z — l)aRT, where aRT represents the bond energy between two monomers [40]. 
For such a model, a reaction is neutral from the point of view of chemical potentials, i.e. when fin+m = Mn + 
when the polymer chain is sufficiently long so that ~ —laRT. From Eq. (13), it follows that 

jL-r ("leq 

RT\n = aRT, (37) 

•^nm 

which then implies a relation between the parameter A defined earlier as the ratio of the rate constants (assumed 
constant) and the parameter a, namely A = 

Since we introduced aRT as bond energy between two consecutive monomers, the simplest choice is to assume 
that leads to a molar enthalphy = —{I — l)aRT, and a molar entropy which is assumed to be negligible 
with respect to the enthalpy part due to bond formation. With this choice, one finds the following contribution of the 
polymer to the enthalpy: 


Hi = -aRT^yi{l - 1)= -aRT{M - cit)). (38) 

;>i 

As discussed previously, at the two fluid level, this should be complemented by the contribution of the solvent to obtain 
the enthalpy H. Similarly, the system entropy S, which contains both contributions, is 

S = -R^yi\nyi. (39) 

/>o 



FIG. 3. System entropy S for the two-fluid version of the String model as a function of time. The various curves represent different 
initial conditions as in Fig. 2. The solvent concentration is co = 1000 and the transition rates are the same as in Fig. 1. 

In Fig. 3, we show this entropy function as a function of time for different initial conditions in units of R. Alt = 0, 
it equals approximately (1 + lnco)/co, and it converges towards the equilibrium value of the entropy at large times. 
For M > 1, the entropy increases monotonically whereas we see that it decreases for M = 1. As in the case of the 
one fluid model, this decrease is not inconsistent since the entropy is not the Fyapunov function. We note that the 
non-monotonicity that was present in the one-fluid model in Fig. 2 for M > 1 is absent in the two-fluid case. 

If the monomers were to assemble in the polymer in a different way, for instance in the form of disks instead of 
linear chains, the standard chemical potentials would be different. In such a case, under similar assumptions as above, 
these chemical potentials would be of the form /i° = — (/ — '/l)aRT, where a is again some constant characteristic 


9 


















of the monomer-monomer and monomer-solvent interaction [40]. The term in a/I represents the contribution of the 
surface energy of the cluster of size 1. This term necessarily implies that the rate constants k~^ must depend 
on n and m in order to satisfy the local detailed balance condition Eq. 13. For such a case, the above derivation of a 
simple exponential for the equilibrium distribution would no longer hold and both the equilibrium and the dynamics 
will would be more complex. 


IV. REVERSIBLE POLYMERIZATION MODEL WITH TWO CONSERVATION LAWS 


As done previously in section II for a reversible polymerization model with only one conservation law, namely 
the total monomer concentration M, we now carry out a similar analysis for a different class of models with two 
conservation laws, namely M and the total concentration of polymers or clusters, c. Clearly, the latter quantity varies 
in time in the String model because some exchange process in Eq. (1) produce clusters of zero length for some n or 
m. In order to construct a model, which conserves the total number of clusters, one needs to forbid such transitions. 
One simple way to achieve this is to consider the kinetic model 

[n] + [m] ->• [n + l] + [m- 1], (40) 

knm 


with the condition n > 1 and m > 2, where the latter inequality precisely prevents the forbidden transitions. It is 
easy to check that now the total concentration of the polymers, c = as well as the total monomer concentration 

M = remain constant in time. Another important observation is that this model is fully reversible even if we 

do not indicate backward reactions explicitly. Indeed, it would be redundant to do so, since backward reactions are 
already included in the forward reactions via an appropriate choice of the indexes (n, m). As done in the previous 
section, we first present a general proof of convergence to equilibrium and then we make contact with thermodynamics 
by introducing chemical potentials in dilute solutions. 

The equation for the rate of change of concentration for polymer size distribution is 

OO OO 

Q — 0(/ 2) ^ '^ [kl — i^n+lCn+lC-l — l knlCnCl\ -f ^ ^ (Tl) 

n—1 m—2 


where the Heaviside function 0(/ — 2) equals 1 for I > 2, and is zero otherwise. 
Assuming again elementary reactions, the entropy production rate E is 


E 


n>l.m>2 


U^nmCnC-m ^m —l,n+lCn+lC77i—l] In 


^nm^n 


1 ,n+1 1 Ctti — 1 


> 0 , 


(42) 


which vanishes when the detailed balance condition holds, Le. in equilibrium; 


U — L 

’^TiTn^ri 


eq eq 


(43) 


Following a procedure similar to that of section II, one can show that, since c is constant, the relative entropy between 
the distribution xi and 


L = R Cl In = Rc Xi In 


xi 

eq 5 


(44) 


is a Lyapunov function. Indeed, this quantity is convex, non-negative (by the inequality In a; < a; — 1), and a mono- 
tonically decreasing function vanishing at equilibrium. This latter property follows from 


dL 

dt 


R Cl In 

i 



(45) 


which using Eq. (41) and the detailed balance condition of Eq. (43) gives 


dL 

dt 


R ^ knlCnClln 
n>l,l>2 


kl — 1 ,n+1 ^71+1Q — 1 

knl^nC-l 


(46) 
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After symmetrizing this sum, we recover that T, = —dL/dt > 0. This result is equivalent to Eq. (6) in presence of the 
additional conservation law c = 0. This system will therefore relax to a unique equilibrium state, where S vanishes. 

We now turn to the two fluid version of the model. As in section IIB, for the two fluids model the molar fraction 
of a polymer of length I is yi{t) = ci{t)IC, where C = c + cq is again a constant. The change in chemical potential 
during the reaction (40) is given by 


Ap — “t” y^m—l Mn 

= /i°+i + Mm-l - - Mm + RT In 

Since at equilibrium Afx = 0, using (43), we get that 


yn-\-iym — l 
ynym 


I, 

0 I 0 0 0 1D'T~'^ 2 / 771 —1 1 H171 

Mn Mm /^n+1 Mm—1 RT In gg RT In — 

2 /n ym ^m — 1 , 


71+1 


The enthalpy change (9) can be written as 


dt 


I I 

^ ^ (/Cm—1.77+1^71+lCm —1 H” 

n>l,m>2 




7>l,m>2 


cy ^ ^ iJ ^ nm ^ nC-m ^m —1,71+1^77+10771 — 1) “h 


7i>l,m>2 


^0 X 


(47) 


(48) 


(49) 


and the entropy change (10) as 

^ = ~ - R\nci) 


l>0 


l>0 


E 

i>l,m>2 


// 7 ^\/0 lO OOiDi ynym \ 

\^nm^n^m ^m—l^n-\-l^n+lCm—l) I ^m —1 ^ 7 i+l ^m -itIn 1 

V [yn+iym-lj/ 


(50) 


Since the entropy production can be rewritten as 


TYu — RT ^ ^ [^i^TimOTiCTTx AJtti—I^ti-I-iCti-I-iCtti —i] In 


^nm^n^m 


i>l,m>2 


2^m —1,71+1 Cti+I Cm— 1 


^ ^ [^i^nmOTiCTTx /u 777 _ 1 ^ 77 _|_lC 77 -(_lC 777 —l] 


i>l,m>2 


M° + Mm - M+1 - Mm-1 + RT^n 


ynVr, 


yn+iym—1 


„ {knmCnCm fcm-l,n-|-lCn-|-l Cm-1) (M 77 H” Mn+l Mm —l); 


(51) 


n>l,m>2 

we can express it, as in (20), as 


/>0 


(52) 


To summarize, we recover exactly the same results as in section II, provided we treat the total polymer concentration 
c as a constant. 
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V. APPLICATION TO THE KINETICS OF GLUCANOTRANSFERASES DPEI AND DPE2 


In this section, we consider the polymerization of glycans by two enzymes studied by Kartal et al. [18], namely 
the glucanotransferases DPEI and DPE2. We show how to construct dynamical models that are compatible with the 
equilibrium polymer length distributions that they found and we study the dynamics of the Shannon entropy for various 
initial conditions. 


A. Kinetics of glucanotransferases DPEI 

Let us assume that the initial condition is not purely made of monomers, since the solution of Ref. [18] becomes 
singular in that limit (see Eq. (4) on P3 of Ref. [18] when the parameter DPin = 1 for instance), and let us construct 
an appropriate dynamics, choosing for simplicity constant rates knm = ^ independent of n and m. 

Using Eq. (41), we have for I > 2, 

dci 

— = k[c (c/+i + Q _1 - 2ci) - CiCi-i + Ciq], (53) 

where the second term forbids transitions from [1] + [Z — 1] to [0] + [/], while the last term forbids transitions from 
[1] + [1] to [0] + [/ + !]. Similarly, for ci, the evolution is 

^ = k[c(c2 - Cl)+ Ci]. (54) 

It is straightforward to verify that this dynamics has two conservation laws, namely J2i>i c; = c and X]z>i = M. 

Introducing the generating function x{z, t) = J2i>i ci{t)z^ as in section III leads once again to a set of equations 
for the dynamics which unfortunately can not be solved analytically. However, it enables us to hnd an explicit solution 
for the equilibrium state: 


= x{z,t ^ oo) 


M — Mz + cz^ 


(55) 


which means that the size distribution ci tends towards the following equilibrium distribution for Z > 1: 


eq _ ^ 

“ M 



(56) 


where the total polymer concentration c is now hxed by the initial condition. We note that the form of the equilibrium 
distribution is the same as that of the String model, but in Eq. (34) is different from c{t = 0), whereas in the 
present DPEI model they are the same. Our equilibrium solution (56) also matches that of Kartal et al. found for the 
polymerization of glycans by glucanotransferases DPEI [18]. In this reference, the authors use the polymer fractions 
xi, where Z stands for the number of linkages in one cluster, rather than our cluster distribution cz. They are related by 
XI = cz+i/c. Their conservation laws therefore read J2i>o ^ J2i>o = DPini — 1, where DPini stands 
for the initial degree of polymerization. The latter is related to c by DPini = M/c, and the relation 1 — c/M = e~^ 
matches Eq. (4) of Kartal et al. 

The Shannon entropy, Eq. (8), at equilibrium and in R units, reads 


„ / X M 

Ssh{t oo) = —In- 

c c 



(57) 


and has the standard form of a mixing entropy. Note that the case of monomer-only initial condition, namely c = 
M = 1, is singular since no evolution is possible from this initial condition according to the present dynamics. In this 
case, the Shannon entropy stays at zero, for all times t, whereas for other initial conditions it increases monotonically 
as shown in Eig. 4. 
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FIG. 4. Left (resp. right): Evolution of the Shannon entropy Ssh for the one fluid DPEl (resp. DPE2) model with time for different 
initial conditions. For both plots, the curves for M — 1,2,3 and 4 correspond to initial conditions of the form ci (0) = 5im- For the 
DPE2 plot on the right, M = 4.6 corresponds to an initial condition of the form c; (0) = OASn + 0.6Sir. Note that when M = 1 
for DPEl and M — 1,2 and 3 for DPE2, the Shannon entropy is zero since no evolution is possible from the initial condition. We 
have chosen the constant rate k = 1. 


B. Application to the kinetics of glucanotransferase DPE2 

As discussed by Kartal et al. [18], the enzyme DPE2 introduces an additional constraint with respect to the enzyme 
DPEl. This additional constraint, which imposes a conservation of the total number of monomers and dimers, reads 
in our notations Ci + C 2 = pc, where p depends on the initial conserved total number of molecules of maltose 
(corresponding to C 2 ) and glucose (corresponding to ci). This additional constraint requires a modification of the 
dynamical evolution equations. We propose the following modification: 

^ = k[c(c 2 - Cl) +cl-cl + C 1 C 3 ], 

dc2 dci 

dt dt 

= k[c(c 4 - C 3 ) - C 2 C 4 + C 1 C 3 + Cg], 

and for I > 4, 

dci 

— = k[c (c /-|-1 + Cl-i — 2ci) — CiC/_i — C 2 Ci+i 

-C3CZ_1 + (ci + C 2 + C3)ci], (61) 

As in the case of DPEl, one can solve the stationary state of this equation by means of generating functions. One 
obtains the following stationary generating function; 


(58) 

(59) 

(60) 


00) 


[ciC2 — CCl — (C1C2 — CC2 + CC3)] Z 

c(z — 1) — ClZ + C2 — C3Z 


(62) 


One obtains from this ci = pc/{l + /) and C 2 = cif = fpcl{l + /) with / = (c — Ci — C 3 )/(c — C 2 ), and for I > 3, 
Cl = £ 3 /*“^ with C 3 = (1 — p)(l — f)c. In other words, for DPE2, the equilibrium distribution is again exponential 
but only for length I > 3, for I < 3 the ratio of C2/C1 for instance does not match the ratio ci+ijci for I > 3. The 
quantity / can be written in terms of p and c only 

M-2c{l-^p) =pcj^+ c{1-p)y^, (63) 

which matches Eq. (S57) obtained by Kartal et al. [18]. Therefore, the equilibrium state (62) is the same as that 
discussed in this reference. 
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We have thus proposed dynamical models reproducing the equilibrium distribution of glycans in presence of DPEl 
or DPE2. The difference between both situations is that DPEl has two conservation laws, namely that of M and of 
c, while DPE2 has a third one corresponding to that of p. As a result, there are more initial conditions of the type 
ci{0) = SiM from which no evolution is possible in DPE2 (M < 3) as compared to DPEl. When this happens, 
Ssh = 0 as shown in Eig. 4. While this forbids initial conditions of pure dimers for instance, no such constraint exists 
for mixtures. Eor instance, an initial mixture of 40:60 of maltose and maltoheptaose considered in [18], corresponding 
to cj(0) = OASii + 0.6(5(7, has p = 0.4 and M = 4.6, and evolves according to DPE2 dynamics as shown in Eig. 4, 
while an initial solution of pure maltose would not. 

The limiting value of the Shannon entropy at long times can be obtained analytically as a function of / and p for 
any initial conditions, but the expression is lengthy and will not be given here. We have checked that it reproduces the 
correct value of the plateaux in Eig. 4. 


VI. CONCLUSION 

In this paper, we have considered two classic models for reversible polymerization in closed systems following the 
mass-action law, one preserving the total polymer concentration and the other one not. In both cases, the entropy 
production can be written as the time derivative of a Lyapunov function which guarantees the relaxation of any initial 
condition to a unique equilibrium satisfying detailed balance. As such, these models could also describe non-chemical 
systems undergoing an aggregation-fragmentation dynamics. 

When considering the polymerization dynamics in dilute solutions, we have shown that a consistent nonequilibrium 
thermodynamics can be established for both models. We find that entropy production is minus the time derivative of 
the nonequilibrium free energy of the system, which is a Lyapunov function and takes the form of a Kullback-Leibler 
divergence between the nonequilibrium and the equilibrium distribution of polymer length. A related result was found 
for the cyclical work performed by chemical machines feeding on polymers in Ref. [41]. Similar relations expressing 
the entropy production as a Kullback-Leibler divergence between the nonequilibrium and equilibrium distributions 
have also been found or used in many studies on Stochastic thermodynamics [42^4]. 

As an application of reversible polymerization models which do not preserve the total polymer concentration, we 
have studied the String model. In this model, the rates of aggregation and fragmentation are constants, which leads to 
an exponential equilibrium distribution of polymer length. At the one-fluid level, we have observed that the Shannon 
entropy is non-monotonic, which is allowed since it differs from the Lyapunov function. At the two-fluid level where 
there is a proper nonequilibrium thermodynamics, no such non-monotonicity arises. 

As an application of reversible polymerization models preserving the total polymer concentration in addition to the 
total number of monomers, we have studied two specific examples named DPEl or DPE2 after Ref. [18]. We have 
shown how to construct dynamics which converge at long times to the expected form and we have discussed the time 
evolution of the Shannon entropy at the one-fluid level. In all the cases, we have been able to And the form of the 
stationary distribution, by applying the method of generating functions. This method is general and also applicable to 
situations where the stationary distribution is a nonequilibrium one [45]. 

Key assumptions of our approach are that we disregarded fluctuations, assumed homogeneous and ideal solutions, 
considered closed systems, and we treated the polymerization reactions as elementary. Each of these assumptions 
could in principle be released and the resulting implications analyzed. Another interesting future direction concerns 
the study of nonequilibrium thermodynamic devices or strategies which can be used to engineer a particular poly¬ 
mer distribution (for instance a monodisperse one) starting from an initial polydisperse one (an exponential one for 
instance). 
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